from EOS_collection import CUBIC_EOS
import numpy as np

#输入对应态参数
Tc = 365
Pc = 4.62e6
w = 0.148
Z = 0.285
R = 8.314
M = 18
Vc = Z * R * Tc / Pc

#输入实际的条件
T = 322
PS = Pc * 10 ** ( 7 * ( 1 + w ) / 3 * ( 1 - Tc / T ) )
phiSL = 1.1
phiSV = 0.1
Tr = T / Tc

#求PR方程的参数
PR = CUBIC_EOS(Pc, Tc, Z, w)
a = PR.PR_a(T)
b = PR.PR_b()

while abs(phiSV - phiSL) > 1e-5:
    #计算气相的体积根
    VSV = R * T / PS
    VSVnew = PR.PR_V(PS, VSV, T, a, b)
    while np.sum ( np.abs(VSVnew - VSV) ) / np.sum ( VSV ) > 1e-5:
        VSV = VSVnew
        VSVnew = PR.PR_V(PS, VSV, T, a, b)

    #计算液相的体积根
    VSL = 2 * b
    VSLnew = PR.PR_V(PS, VSL, T, a, b)
    while np.sum ( np.abs(VSLnew - VSL) ) / np.sum ( VSL ) > 1e-5:
       VSL = VSLnew
       VSLnew = PR.PR_V(PS, VSL, T, a, b)

    #计算压缩因子
    ZSV = PS * VSV / R / T
    ZSL = PS * VSL / R / T

    #计算气相的逸度系数
    phiSV = PR.PR_PHI(PS, VSV, T, a, b, ZSV)

    #计算液相的逸度系数
    phiSL = PR.PR_PHI(PS, VSL, T, a, b, ZSL)

    PS = PS * ( 1 - (phiSV - phiSL) / (ZSV - ZSL))

print(f"饱和蒸汽压为{PS*1e-3} kPa")